Map vulnerability and impacts
Compiled on Tue Jun 22 17:56:49 2021 by ohara
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
library(raster)
library(oharac)
oharac::setup()
source(here('fb_slb_fxns.R'))1 Summary
Combine vulnerability scores with AquaMaps species information and functional groups.
2 Data
See individual trait scripts.
3 Methods
3.1 Gather and assemble data
spp_fxn_groups <- read_csv(here('_data/functional_groups_from_traits.csv'))
spp_info <- get_am_spp_info()
spp_by_gp <- spp_info %>%
dplyr::select(am_sid, sciname) %>%
inner_join(spp_fxn_groups, by = c('sciname' = 'species')) ### 18285 spp matched in fxn gps
spp_vuln <- get_spp_vuln() %>%
select(sciname = species, stressor, score, sd_score) %>%
distinct()
spp_vuln_am <- spp_info %>%
inner_join(spp_vuln, by = c('sciname'))
nspp <- spp_vuln_am$sciname %>% n_distinct()
### 10016 species with vulnerability scores in the AquaMaps dataset overall
### (will improve when AquaMaps names are rectified with WoRMS)
spp_vuln_by_gp <- spp_by_gp %>%
inner_join(spp_vuln, by = c('sciname'))
nspp <- spp_vuln_by_gp$sciname %>% n_distinct()
### 7985 species with both functional groups and vulnerability scores. Of
### these, the highest group-gapfill level is ~3.2, not so bad...3.2 Examine vulnerability by stressor and functional group
if(file.exists(here('figs/vuln_plot_by_group.png'))) {
message('plot exists, not going to regenerate it!')
} else {
spp_vuln_variation <- spp_vuln_by_gp %>%
group_by(stressor) %>%
mutate(diff = score - mean(score)) %>%
group_by(p_group) %>%
mutate(gp_lbl = paste0(p_group, ' (n = ', n_distinct(sciname), ')')) %>%
ungroup()
spp_vuln_var_summary <- spp_vuln_variation %>%
group_by(stressor, gp_lbl) %>%
summarize(mean = mean(diff),
sd = sd(diff)) %>%
ungroup()
x <- ggplot() +
theme_minimal() +
theme(axis.text = element_text(size = 7)) +
geom_jitter(data = spp_vuln_variation,
aes(x = stressor, y = diff), color = 'grey70', alpha = .1) +
geom_hline(yintercept = 0, color = 'red', size = 1) +
geom_linerange(data = spp_vuln_var_summary,
aes(x = stressor, ymin = mean - sd, ymax = mean + sd),
color = 'black', size = .5) +
geom_point(data = spp_vuln_var_summary, shape = 21, size = 1.5,
aes(x = stressor, y = mean), color = 'black', fill = 'yellow') +
geom_hline(yintercept = 0, color = 'red', size = .25) +
coord_flip() +
theme(axis.title.y = element_blank()) +
facet_wrap(~gp_lbl) +
labs(y = 'difference (vuln score - mean(vuln score))')
ggsave(here('figs/vuln_plot_by_group.png'), height = 8, width = 8)
}
3.3 Map CHI for one species
Select a species (or a small handful) and calculate cumulative impacts across range.
- Dermochelys coriacea (leatherback turtle): Rep-3437
- Poroderma pantherinum (Leopard catshark): Fis-23221
- Megaptera novaeangliae (humpback whale): ITS-Mam-180530
- Thunnus orientalis (Pacific bluefin tuna): Fis-123736
- Clupea harengus (Atlantic herring): Fis-29344
While we have stressor data from BD/CHI at 10km x 10km resolution, species distributions are at 0.5° resolution… for now let’s use some test rasters (from setup_stressors.Rmd) of the stressors in their most recent available year.
3.3.1 The stressors
str_maps <- data.frame(f = list.files(here('_spatial/stressors_hcaf'), full.names = TRUE)) %>%
mutate(str_chi = str_remove(basename(f), '_[0-9]{4}.+'),
year = str_extract(basename(f), '[0-9]{4}') %>% as.integer()) %>%
full_join(read_csv(here('_raw/strs_vuln_to_chi.csv')) %>%
mutate(str_chi = str_split(str_chi, ';')) %>%
unnest(str_chi)) %>%
filter(!is.na(str_vuln) & !is.na(str_chi) & str_vuln != 'biomass_removal')
str_stack <- raster::stack(str_maps$f) %>%
setNames(str_maps$str_vuln)
### fill NAs with zero for math reasons
values(str_stack)[is.na(values(str_stack))] <- 0
plot(str_stack)
spp_vec <- c('Dermochelys coriacea (leatherback turtle)' = 'Rep-3437',
'Carcharodon carcharias (great white shark)' = 'Fis-23071',
'Megaptera novaeangliae (humpback whale)' = 'ITS-Mam-180530',
'Thunnus orientalis (Pacific bluefin tuna)' = 'Fis-123736',
'Clupea harengus (Atlantic herring)' = 'Fis-29344')
test_vuln_df <- spp_vuln_by_gp %>%
filter(am_sid %in% spp_vec) %>%
select(am_sid, sciname, p_group, stressor, score, sd_score) %>%
filter(stressor %in% str_maps$str_vuln)
test_maps_df <- get_am_spp_cells(prob = 0) %>%
filter(am_sid %in% spp_vec) %>%
mutate(presence = 1)
for(i in seq_along(spp_vec)) {
### i <- 2
spp <- spp_vec[i]
message('processing impacts for ', names(spp_vec)[i])
cat(sprintf('\n\n### %s \n', names(spp_vec)[i]))
tmp_spp_range_rast <- test_maps_df %>%
filter(am_sid == spp) %>%
map_to_hcaf(which = 'presence')
tmp_spp_vuln_df <- test_vuln_df %>%
filter(am_sid == spp) %>%
filter(score > 0)
impact_list <- vector('list', length = nrow(tmp_spp_vuln_df)) %>%
setNames(tmp_spp_vuln_df$stressor)
for(j in 1:nrow(tmp_spp_vuln_df)) {
### j <- 1
s <- tmp_spp_vuln_df$stressor[j]
v <- tmp_spp_vuln_df$score[j]
v_rast <- tmp_spp_range_rast * v
s_rast <- str_stack[[s]]
i_rast <- v_rast * s_rast
impact_list[[j]] <- i_rast
}
impact_stack <- stack(impact_list)
impact_stack[['chi']] <- calc(impact_stack, sum)
plot(impact_stack)
}3.3.2 Dermochelys coriacea (leatherback turtle)
3.3.3 Carcharodon carcharias (great white shark)
3.3.4 Megaptera novaeangliae (humpback whale)
3.3.5 Thunnus orientalis (Pacific bluefin tuna)
3.3.6 Clupea harengus (Atlantic herring)
spp_vec <- spp_vuln_am %>%
filter(stressor %in% str_maps$str_vuln) %>%
.$am_sid %>%
unique()
spp_vuln_df <- spp_vuln_am %>%
filter(am_sid %in% spp_vec) %>%
select(am_sid, sciname, stressor, score, sd_score) %>%
filter(stressor %in% str_maps$str_vuln)
spp_maps_df <- get_am_spp_cells(prob_cut = 0) %>%
filter(am_sid %in% spp_vec) %>%
mutate(presence = 1)spp_chi_fstem <- here_anx('spp_impacts_hcaf/%s_chi_hcaf.csv')
### unlink(here_anx('spp_impacts_hcaf/*.*'))
n_done <- list.files(here_anx('spp_impacts_hcaf'), pattern = 'chi_hcaf.csv') %>%
n_distinct()
tmp <- parallel::mclapply(spp_vec, mc.cores = 12,
FUN = function(spp) {
spp_chi_file <- sprintf(spp_chi_fstem, spp)
if(file.exists(spp_chi_file)) {
# message('File ', basename(spp_chi_file), ' exists... skipping!')
} else {
message('Processing impacts for ', spp)
tmp_maps_df <- spp_maps_df %>%
filter(am_sid == spp)
tmp_spp_range_rast <- tmp_maps_df %>%
map_to_hcaf(which = 'presence')
tmp_spp_vuln_df <- spp_vuln_df %>%
filter(am_sid == spp) %>%
filter(score > 0)
impact_list <- vector('list', length = nrow(tmp_spp_vuln_df)) %>%
setNames(tmp_spp_vuln_df$stressor)
for(j in 1:nrow(tmp_spp_vuln_df)) {
### j <- 1
s <- tmp_spp_vuln_df$stressor[j]
v <- tmp_spp_vuln_df$score[j]
v_rast <- tmp_spp_range_rast * v
s_rast <- str_stack[[s]]
i_rast <- v_rast * s_rast
impact_list[[j]] <- i_rast
}
impact_stack <- stack(impact_list)
impact_stack[['chi']] <- calc(impact_stack, sum)
impact_df <- as.data.frame(values(impact_stack)) %>%
mutate(across(everything(), ~round(.x, 5))) %>%
mutate(loiczid = 1:n()) %>%
filter(!is.na(chi)) %>%
### this works because all "present" cells have a numeric value (incl 0)
left_join(tmp_maps_df, by = 'loiczid') %>%
select(-am_sid, -presence)
write_csv(impact_df, spp_chi_file)
}
})3.4 Calculate ocean-area-weighted mean CHI for each species…
a_df <- get_hcaf_info() %>%
select(loiczid, ocean_area) %>%
distinct()
spp_chi_fs <- list.files(here_anx('spp_impacts_hcaf'), pattern = 'chi_hcaf.csv', full.names = TRUE)
spp_ids <- str_remove(basename(spp_chi_fs), '_chi_hcaf.csv')
mean_chi_list <- parallel::mclapply(spp_chi_fs, mc.cores = 16,
FUN = function(f) {
### f <- spp_chi_fs[1]
spp <- str_remove(basename(f), '_chi_hcaf.csv')
spp_chi <- data.table::fread(f) %>%
mutate(am_sid = spp) %>%
left_join(a_df, by = 'loiczid') %>%
group_by(am_sid) %>%
summarize(area_km2 = sum(ocean_area),
mean_chi = sum(chi * ocean_area) / area_km2)
if(nrow(spp_chi) == 0) {
spp_chi <- data.frame(area_km2 = NA, mean_chi = NA)
}
}) %>%
setNames(spp_ids)
mean_chi_df <- bind_rows(mean_chi_list, .id = 'am_sid')
write_csv(mean_chi_df, here('int/mean_chi_check.csv'))mean_chi_df <- read_csv(here('int/mean_chi_check.csv')) %>%
left_join(spp_vuln_by_gp %>% select(am_sid, p_group) %>% distinct(), by = 'am_sid')
ggplot(mean_chi_df, aes(x = log10(area_km2), y = mean_chi, color = p_group)) +
geom_point(show.legend = FALSE) +
geom_hline(yintercept = 0) +
facet_wrap( ~ p_group) +
theme_minimal()